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We study numerically phase separation in a binary fluid subject to an applied shear flow in two 
dimensions, with full hydrodynamics. To do so, we introduce a mixed finite-differencing/spectral 
simulation technique, with a transformation to render trivial the implementation of Lees-Edwards 
sheared periodic boundary conditions. For systems with inertia, we reproduce the nonequilibrium 
steady states reported in a recent lattice Boltzmann study. The domain coarsening that would 
occur in zero shear is arrested by the applied shear flow, which restores a finite domain size set by 
the inverse shear rate. For inertialess systems, in contrast, we find no evidence of nonequilibrium 
steady states free of finite size effects: coarsening persists indefinitely until the typical domain size 
attains the system size, as in zero shear. We present an analytical argument that supports this 
observation, and that furthermore provides a possible explanation for a hitherto puzzling property 
of the nonequilibrium steady states with inertia. 

PACS numbers: 64.75. +g, 47.11. Qr. 



I. INTRODUCTION 

When an initially homogeneous mixture of two incom- 
pressible fluids (A and B) undergoes a deep temperature 
quench into the spinodal regime, it becomes unstable 
with respect to spatial fluctuations in the composition 
fleld (l){x,t). The mixture then phase separates into well 
deflned domains of A-rich and B-rich fluid, which, after 
a rapid initial transient, attain local equilibrium within 
each domain. The system remains globally out of equilib- 
rium on long timescales, however, due to the excess en- 
ergy that resides in the interfaces between the domains. 
Here we consider the maximally symmetric case of a 50:50 
mixture of two mutually phobic fluids with matched vis- 
cosities and densities. 

Following the initial transient of their formation, the 
domains slowly coarsen in time through the action of the 
surface tension in the interfaces that separate them. (For 
reviews of phase ordering kinetics, see Refs. In 
this way, the excess interfacial energy of the system pro- 
gressively relaxes towards its minimal equilibrium value. 
This coarsening process proceeds through three distinct 
regimes that are successively dominated by diffusive, vis- 
cous and inertial dynamics. In the limit of an infinite 
system size A — s- cx) (taken first), the typical domain size 
perpetually increases without bound: the system never 
globally equilibrates, even in the limit of infinite time 
t — > oo (taken second). For any finite system size, coars- 
ening is in practice eventually cutoff once the typical do- 
main size L(t) attains the system size A. 

Besides these systems that remain out of equilibrium 
as they slowly relax towards a fully phase separated 
state, another class of systems comprises those that are 
continuously driven out of equilibrium by the external 
application of a steady shear flow. In this paper, we 
consider systems that are both undergoing phase separa- 
tion and simultaneously subject to an applied shear flow. 
They thus combine both of the nonequilbrium features 



just described. The main question that we address is 
whether shear interrupts domain coarsening and restores 
a nonequilibrium steady state with a typical domain size 
L(7^^) set by the inverse of the applied shear rate 7; or 
whether coarsening persists indefinitely, up to the system 
size, as in zero shear. 

Despite previous experimental [4l [5l [6l |7l [8] [9] [TO] 
[HI Ea Ha [H], numerical [15l [TBI [17l UHl [HI Efll EH 
[lll[131lllll51llSl[lZl[Ml[Ml[3nj[3I] and theoreti- 
cal [TSl[Tnj|inillIl[2SJ[33[Ml[Ml[S51[MJ[S3 work, 

this question remained open until the recent simulation 
studies of Stansell et al. in Refs [351 Using Lattice 
Boltzmann techniques, they gave convincing evidence for 
the formation of nonequilibrium steady states, with do- 
mains of a finite size set by the inverse shear rate, in- 
dependent of the system size. The domain morphology 
inherits the anisotropy of the applied shear flow, and is 
therefore characterised by two lengthscales Ly and L±, 
respectively describing the major and minor domain axes. 
A remarkable achievement of Refs. [551 132] was the ex- 
ploration, by judicious parameter steering, of a range of 
inverse shear rates spanning six decades on a scaling plot. 
For reasons discussed below this range would, a priori, be 
expected to encompass two different regimes, separately 
dominated by viscous and inertial dynamics. Surprising, 
therefore, is the observation of an apparently single scal- 
ing with shear rate for each of L|| and L±, across all six 
decades. 

All the simulations reported in Refs. [38l[39] have non 
zero fluid inertia. Indeed, even when attempting to ac- 
cess the limit of zero inertia, a small but non-zero inertia 
remains a practical requirement of the Lattice Boltzmann 
technique. In this paper, we investigate phase separation 
under shear in systems that are strictly inertialess. To do 
so, we apply a different simulation technique to this prob- 
lem, comprising finite differencing combined with Fourier 
spectral methods. We also use a convenient transforma- 
tion to render trivial the implementation of sheared pe- 
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FIG. 1: Flow geometries in the (a) unsheared and (b) sheared cases. The crossed connected by thin dotted hnes show 
representative pairs of points connected by periodic boundary conditions. 



riodic boundary conditions I40j. 

Our main numerical result will be that, in truly incr- 
tialess systems, coarsening persists indefinitely, up to the 
system size, despite the external shear flow. We will also 
present an analytical argument that supports this obser- 
vation. A simple extension to the same argument pro- 
vides a possible explanation for the existence of a single 
scaling for each of Ly , L_l across all six decades of scaled 
shear rate in the simulations of Ref . [35| , with inertia. 

We start by introducing the model equations, flow ge- 
ometry and choice of units in Sees. |ll III and IV re- 
spectively. In Sec. |V] we briefly outline our numerical 
method, of which further details are given in the Ap- 
pendices. The length and timescales that characterise 



demixing are presented in Sec. VI leading to a discus- 
sion in Sec. |VII| of the choice of parameter values in our 
simulations. In Sec. |VIII| we present our numerical re- 
sults, starting in Sec. VIII A with the case of coarsening 
in zero shear, before considering an applied shear flow 
with and without inertia in Sees. IVIIIBI and I VIII D I re- 
spectively. Se c. |VIII C| contains a linking discussion. The 
results in Sec. VIII D are (we believe) novel. In contrast, 
our aim in Sees. VIII A| and | VIII 5] is simply to reproduce 
the behaviour seen in previous lattice Boltzmann stud- 
ies [351 mi 112 > thereby gaining confidence in our own 
simulation method. In Sec. IX we present an analytical 
argument that supports our numerical observations, as 
well as those of Ref. [3S]. Sec. [X] contains a summary 
and an outlook to future work. 



II. GOVERNING EQUATIONS 

The fluid velocity fleld v{x,t) and pressure p{x,t) are 
governed by the Navier-Stokes equation 

p{dtv + v - Vv) = 7jV'^v- <j)Vfi-Vp (1) 

together with the incompressibility constraint 

V ■v = 0. (2) 



Here p is the fluid density and 77 the viscosity. p{x,t) 
is the thermodynamic part of the pressure tensor and 
fi{x,t) is a chemical potential, deflned below. In what 
follows, we shall consider two dimensional (2D) flow in 
the X — y plane with velocity components Vx,Vy. 

To eliminate the pressure, we take the curl of Eqn. [l] 
The z component of the resulting equation is 

p{dtLU + v-\/u;)^-qV'^u;~[V A(t)Vp]-z (3) 

in which the vorticity to obeys 

V Av = ujz. (4) 

To ensure that the incompressibility constraint Q is au- 
tomatically obeyed, we define a streamfunction ijj via 



v V A tpz, 
which is related to the vorticity as follows: 



(5) 



(6) 



The dynamics of the compositional order parameter 
(f>{x, t) are prescribed by an advection-diffusion equation 
of Cahn-HiUiard type (see Refs. 0113]) 



dt<j) + v -S/cj)^ M'^^p. 



(7) 



Here M is the mobility, assumed constant, controlling the 
rate of intermolecular diffusion. The chemical potential 



p = G 



1) - kV^ 



(8) 



in which G is a positive constant with the dimensions 
of stress, k is also a positive constant, controlling the 
characteristic width I of the interface between domains: 

The surface tension of the interface is given by 

.^^^ (10) 
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III. FLOW GEOMETRY 



NUMERICAL METHOD 



In the absence of shear we consider a square domain 
of size Ax — Ay in the x — y plane (Fig. Ilk). We adopt 



(11) 
(12) 



periodic boundary conditions in each direction: 

(l){x^O,y) = (f>{x = A.j„y) V y, 
(l){x,y^O) = (f){x,y = Ay) V x, 



and similarly for Vx,Vy^ij} and uo. 

To study shear, we consider a system that is rectangu- 
lar with Ax = 2 Ay . Conceptually, a shear flow of rate 7 is 
applied by moving the upper boundary to the right at a 
constant speed 7Aj,, as shown in Fig.fTb (left), so that the 
velocity field then has a constant aflme part ^yx and an 
additive fluctuating contribution v{x^ t). In practice, this 
is implemented via Lees-Edwards sliding periodic bound- 
ary conditions (Fig. [TJi, right), such that the system has 
no physical borders but instead comprises an infinitely 
repeating array of images arranged in horizontally slid- 
ing layers. The boundary conditions are then 

4>{x = {),y) = 4>{x^Ax,y) V y, (13) 
0(a;,y==O) = (l){x + -^Ayt.y = Ay) V x. (14) 

These also apply to Vx,Vy, as well as the correspond- 
ing fluctuating contribution to the streamfunction ip and 
vorticity uj. 

IV. CHOICE OF UNITS 

The model equations contain five parameters: 
p, f], M, G and I. (We must specify two of G, I, k and cr, 
and here have chosen G and /.) In addition, we must also 
specify the system's dimensions Ax, Ay and the applied 
shear rate 7, making eight parameters in all. Three of 
these can be eliminated by choosing appropriate units for 
length, time and stress. Accordingly, we measure lengths 
in units of the vertical system size Aj^ = 1; stresses in 
units of G = 1; and times in units of the characteristic 
"microscopic" time for diffusion across an interface, 



Tn = = i. 



(15) 



In these units, the governing equations comprise 
Eqns. |3] [5] and [6] (unchanged) together with 

V -^(1) = fV^fi (16) 



and 



(17) 



This leaves the rescaled density p, viscosity 77, interfacial 
width aspect ratio A^^ and applied shear rate 7 as the 
five control parameters. In our numerical work we fix A^ 
and /, varying p, 77 and 7 between runs. 



In the laboratory frame, the numerical implementa- 
tion of sheared periodic boundary conditions is rather 
involved. We now discuss a convenient transformation 
that renders it trivial |1D]. We do so here in outline only, 
referring the interested reader to appendix [A] for details. 

The transformation comprises two steps. In the first, 
the velocity field is expressed as the sum of a constant 
affine part and a fluctuating contribution: 

v{x,t) — jyx + v{x,t). (18) 

In the second step, we transform to the cosheared frame 

(x, y, t) {x ^x~ ity, y = y, t' = t), (19) 

defining for convenience the cosheared gradient operator 

-^c^-xdx' + y i-itdx' + dy>). (20) 

As shown in Appendix [X] the final transformed equation 
set, dropping the dashes for clarity, is then: 



(21) 



together with 



pdtUj+p{dyipdxU! - dxipdyuj) = r]Vluj-[dx(f>dyp - dycj)dxp] 



dt^ + (dy^dx^ - dx^dy^) = PvIp, 



(22) 
(23) 

(24) 



A priori, the bracketed expressions in Eqns. [22] and [23] 
contain terms in the applied shear rate 7. However in 
each bracket these are equal and opposite, and so cancel. 
An equivalent statement is that the bracketed expressions 
are invariant under shear. 

This transformation considerably simplifies the imple- 
mentation of an applied shear flow. Indeed, looking at 
Eqns. [21] to |24] we see that the only effect of shear is 
to renormalise the gradient operator V — *■ Vc. Further- 
more, the dynamical variables 0, ip and lo are now subject 
to ordinary periodic boundary conditions: 

(l){x = 0,y) = (j){x^Ax,y) V y, 

cj){x,y = 0) = cp{x,y^Ay) V x, (25) 

and similarly for ip and uj. This will allow us to use 
Fourier transforms in our numerical algorithm, as de- 
scribed below. 

Under the transformation described so far, the relative 
shear of the laboratory and cosheared frames becomes 
large at long times, diverging at a constant rate 7. This 
is clearly expected to give rise to numerical instabilities. 
To circumvent this problem, once the relative shear s{t) 
attains Ax/ 2 Ay we perform by hand an instantaneous 
shift s{t) — > s{t) — Ax/ Ay. In this way, the function s{t) 
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is bounded between s = —A.j;/2Ay and s = Ax/2Ay, com- 
prising a sawtooth with sections of constant slope 7 con- 
nected by negative step discontinuities of height Aj./Ay 
at equal time intervals At — Ax/^fAy. Because this shift 
involves moving the upper wall by exactly one multiple 
of the system length, the periodic boundary conditions 
(25 1 are unaffected. For times between these shifts, the 



sole effect of this modification appears in the cosheared 



gradient operator, which, in place of (20 1, is now 



Vc = xdx + y {-s{t)dx + dy). 



(26) 



The special case of zero shear is trivially achieved by 
setting s{t) = OWt. 

In appendix [B] we discuss the numerical algorithm 
used to study the dynamical evolution of (j){x,t),u!{x,t) 
and ip{x, t), as specified by Eqns. 21 to 24 and 26 Read- 



ers who are not interested in these issues can skip straight 
to Sec. IVII without loss of thread. 



VI. LENGTHSCALES AND TIMESCALES 

Following a deep temperature quench of an unsheared 
system into the two phase regime, well defined domains 
of each phase form, separated by sharp interfaces. See 
Fig. [3] for example, in which the white (resp. black) 
patches are domains of A (resp. B)-rich fluid. After a 
non-universal transient associated with their initial for- 
mation, these domains progressively grow in size through 
the action of surface tension, passing through regimes 
that are successively dominated by diffusive, viscous and 
inertial dynamics. During this process of domain coars- 
ening, the dynamical scaling hypothesis [H 33] |35] states 
that any structural lengthscale L{t) characterising the 
typical domain size {e.g., as measured by some moment 
of the structure factor) should depend in the same way 
on time t as any other {e.g., as measured by the total 
amount of interface present in the system) . Within each 
regime, dimensional analysis can be used to construct the 
functional form of L{t) out of the model parameters that 
are relevant to that regime |331 US] . 

• In the diffusive regime, the system is unaware of the 
density p and the viscosity rj. Out of the remaining 
parameters M,a, and the time t we can construct 
only one lengthscale: 

Ld(<) = (Ma^)l/^ (27) 

In our units Ld = alt^/^ with a = (2^2/3)^/^. 

• In the viscous hydrodynamic regime, the system is 
unaware of M and p. From the remaining param- 
eters rj, a and t we have 



Ly{t) = 

In our units Ly = a^lt/rj. 



(28) 



• In the inertial hydrodynamic regime, the system is 
unaware of M and rj. From the remaining param- 
eters p, a and t we have 



Li{t) = 



at' 



In our units L\ = a{l t^ j p) 



1/3 



1/3 



(29) 



By equating the growth laws Ld = ^Vi we can construct 
the characteristic lengthscale Ldv a-t which we expect a 
crossover from diffusive to viscous dynamics. Similarly 
by equating Ly = -^i we find the lengthscale Lyi for 
crossover from viscous to inertial dynamics. Together 
with the "microscopic" interfacial width / and the system 
size Ax, Ay, we then have the following basic Icngthscales 

• The interfacial width /. 

• The characteristic lengthscale for crossover from 
diffusive to viscous hydrodynamic coarsening: 



In our units Ldv ~ 



(30) 



• The characteristic lengthscale for crossover from 
viscous to inertial hydrodynamic coarsening: 



Lyi — 



pa' 



(31) 



In our units Lyi irf / 2^/2 pi. 

• The system size A^;, Aj,. In our units Aj, = 1 always. 
Corresponding to these are the following timescales 

• The microscopic timescale for diffusion across an 
interface 



To = 



2V2 

"3~iWCT 



1. 



(32) 



• The characteristic timescale for crossover from dif- 
fusive to viscous hydrodynamic coarsening 



2%/2 



(33) 



• The characteristic timescale for crossover from vis- 
cous to inertial hydrodynamic coarsening 



VI 



pa^ 



9 

8^' 



(34) 



• The characteristic time at which the domain size 
attains the system size in coarsening 



system ■ 



(35) 
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So far we have discussed the growing structural length- 
scale L(t) in general terms, without any specific defini- 
tion. In fact, there exist many possible measures of the 
characteristic domain size. In zero shear we use the fol- 
lowing one 



/ dq^ J dqyS{q) 
J dq^ J dqy\q\-^S{q) 



(36) 



defined via the structure factor S{q), which, as usual, is 
the 2D Fourier Transform of (l){x,y). 

The discussion so far in this section has related to un- 
sheared systems. Under an applied shear flow the do- 
main morphology becomes anisotropic. See Fig. |6] for 
example. Accordingly, we consider the following ma- 
trix mm 



D. 



I J dx J dy da 
J dx J dy (fy^ 



(37) 



the reciprocal eigenvalues of which give us two length- 
scales L|| , Lj_ characterising the long and short principal 
axes of the domain morphology. 



VII. PARAMETER STEERING 



As discussed in Sec. |IV| the physical control parame- 
ters that must be prescribed in any simulation run are 
the fluid density p and viscosity rj, the interfacial width 
I, the applied shear rate 7, and the system's aspect ra- 
tio Ax- (Recall that we have set Ay = 1, G = 1 and 
To = 1.) We must also specify the number of numerical 
mesh points and Ny, and the numerical timestep St- 
We now discuss appropriate choices of these parameter 
values that will allow us to access the physical regimes of 
interest. 

In the absence of an applied shear flow, 7 — 0.0, our 
aim is to study coarsening of the isotropic domain mor- 
phology. See Fig. |3] for example. In each of these runs, 
we consider a square simulation box with A^; = A^ = 
A = 1, in our units. In any given run, at any time t, our 
aim is to ensure a separation of the four lengthscales 



(5 < ? < L{t) < A. 



(38) 



Here S = 1/N is the mesh size, prescribed by the recip- 
rocal of the number TV = = Ny of numerical mesh 
points in each spatial dimension. Recall that I is the 
width of the interface between domains; L(t) is the grow- 
ing domain size; and A = 1 the system size. We thereby 
restrict ourselves to physical lengthscales I and L{t) that 
lie between the mesh and system size. To make this win- 
dow as large as possible, we take 6 as small as possible, 
using as many grid points iV^ — as is numerically 
feasible. The results presented below have N — 512, 
checked for convergence to the limit ^ 00 at flxcd 
I against N = 1024 (not shown). We set the interfa- 
cial width I to the minimum for which interfaces are still 



fully resolved by this grid, taking I = 0.00156, checked 
for convergence to the limit l/A ^ against I = 0.0025 
and I = 0.005 (not shown). 

In each simulation run, this leaves a window of approx- 
imate size 0.05 < L{t) < 0.25 in which the domain size 
L{t) is much larger than the width I of the interface be- 
tween domains, and much smaller than the system size 
A = 1. Accordingly, we are only able to study a small 
piece of the full curve L(t) in each run. In different runs, 
therefore, we vary the viscosity 77 and density p to en- 
sure that we cover all regimes of interest. As discussed 
in Refs. |l6l|47], we are then able to construct compos- 
ite scaling curves L{t)/L-QY Etnd L(t)/Lvi, each spanning 
several decades of scaled length and time. 

Below we present results for two different series of 
runs. 



In the first, R028u to R032u in table (a) 



we 

explore the viscous and inertial hydrodynamic regimes, 
and the crossover between them, by varying the crossover 
lengthscale iyi- These are in fact the parameter values 
used previously by Kendon and coworkers in their lattice 
Boltzmann simulations f42', converted into our units. 

In the second series of runs, DVlu to DV6u in table |(b)] 
we explore the diffusive and (again) viscous hydrody- 
namic regimes, and the crossover between them. Accord- 
ingly, we set p — Q such that Lyi is infinite and the iner- 
tial regime is never attained. Across the different runs we 
vary the predicted crossover lengthscale Ldv — ^y/v by 
sweeping the viscosity in the range 77=10^^---?7=10^. 
For the largest viscosity values, the system explores the 
diffusive regime. For the smallest values, the system 
passes into the viscous regime as soon as well defined 
domains have formed, without any identifiable regime of 
diffusive domain growth. 

To study the effect of an applied shear flow, we con- 
sider the same sets of parameter values as in zero shear 
but now with non-zero values of 7. In anticipation of 
the anisotropy that shear will induce, we also double the 
length of cell in the flow direction, along with to main- 
tain a constant density of mesh points. See tables (c) 



and (d) as just discussed, the parameters of any sheared 
run (R032s, for example) are the same as those of the cor- 
responding unsheared run (R032u) , apart from the values 
of 7, A3, and N.j;. The parameters of runs R028s, R022s, 
R029s, R020s, R030s, R019s and R032s are those used in 
the lattice Boltzmann study of Ref. [3S], converted into 
our units. The appearance of the new sets R028b, R022b 
and R029b will be explained in Sec. [VITTB] below. 

Assuming for the moment - in some cases incorrectly, 
as we shall show below - that any sheared system will 
attain a steady state with a typical domain size set by 
the reciprocal shear rate, our aim would then be to con- 
struct scaling plots L(l/7) analogous to those of L{t) for 
the zero-shear coarsening regime in Fig. |2] To obtain as 
much information as in Fig. [2] however, we would need to 
run for very many values of the shear rate 7 to produce 
corresponding near continuous segments of L versus I/7. 
This is prohibitive within available computing resources. 
Accordingly, we ran for just one shear rate for each set of 



6 



(a)Zero applied shear. Viscous hydrodynamic to inertial hydrodynamic regime. 
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V 


P 


I 


7 






Ny 


5t 


Ldv 


Lvi 


R028u 


0.111 


136.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000519 


0.0616 


R022u 


0.196 


2510.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000690 


0.0104 


R029u 


0.0467 


893.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000337 


0.00166 


R020u 


0.0785 


16180.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000437 


0.000259 


R030u 


0.0122 


6250.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000172 


0.0000162 


R019u 


0.00876 


32814.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.000146 


0.00000159 


R032u 


0.00391 


20067.0 


0.00156 


0.0 


1.0 


512 


512 


0.008 


0.0000975 


0.000000518 


(b)Zero applied shear. Diffusive to viscous hydrodynamic regime. 
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Ny 


St 




Lvi 


DVlu 


1000.0 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.0494 


oo 


DV2u 


100.0 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.0156 


oo 


DV3u 


10.0 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.00494 


oo 


DV4u 


1.0 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.00156 


oo 


DV5u 


0.3 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.000855 


oo 


DV6u 


0.1 


0.0 


0.00156 


0.0 


1.0 


512 


512 


0.016 


0.000494 


oo 


(c)Applicd shear, with inertia. 
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Ny 


St 




Lvi 


R028s 


0.111 


136.0 


0.00156 


0.0765 


2.0 


1024 


512 


0.008 


0.000519 


0.0616 


R028bs 


0.444 


8704.0 


0.00156 


0.0765 


2.0 


1024 


512 


0.008 


0.00104 


0.0154 


R022s 


0.196 


2510.0 


0.00156 


0.0205 


2.0 


1024 


512 


0.008 


0.000690 


0.0104 


R022bs 


0.339 


22590.0 


0.00156 


0.0355 


2.0 


1024 


512 


0.008 


0.000908 


0.00349 


R029s 


0.0467 


893.0 


0.00156 


0.0341 


2.0 


1024 


512 


0.008 


0.000337 


0.00166 


R029bs 


0.0809 


8037.0 


0.00156 


0.0591 


2.0 


1024 


512 


0.008 


0.000444 


0.000554 


R020s 


0.0785 


16180.0 


0.00156 


0.0256 


2.0 


1024 


512 


0.008 


0.000437 


0.000259 


ROSOs 


0.0122 


6250.0 


0.00156 


0.0410 


2.0 


1024 


512 


0.004 


0.000172 


0.0000162 


R019s 


0.00876 


32814.0 


0.00156 


0.0251 


2.0 


1024 


512 


0.004 


0.000146 


0.00000159 


R032s 


0.00391 


20067.0 


0.00156 


0.051 


2.0 


1024 


512 


0.004 


0.0000975 


0.000000518 


(d)Apphed shear, without inertia. 
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A^ 




Ny 


St 




Lvi 


DVls 


1000.0 


0.0 


0.00156 


0.01 


2.0 


1024 


512 


0.016 


0.0494 


oo 


DV2s 


100.0 


0.0 


0.00156 


0.01 


2.0 


1024 


512 


0.016 


0.0156 


oo 


DVSs 


10.0 


0.0 


0.00156 


0.01 


2.0 


1024 


512 


0.016 


0.00494 


oo 


DV4s 


1.0 


0.0 


0.00156 


0.01 


2.0 


1024 


512 


0.016 


0.00156 


oo 


DVSs 


0.3 


0.0 


0.00156 


0.01 


2.0 


1024 


512 


0.016 


0.000855 


oo 


DV6s 


0.1 


0.0 


0.00156 


0.03 


2.0 


1024 


512 


0.016 


0.000494 


oo 



TABLE I: Parameter values used in simulation runs. 



(other) parameter values, marking the scaled reciprocal 
shear rates recorded in tables |(c)| and |(d)| with vertical 
arrows in Fig. [2] 



VIII. NUMERICAL RESULTS 

We now present the results of our numerical simula- 
tions. We start in Sec. IVIII Al with the case of coars- 



ening in zero shear, before considering an applied shear 
flow with and without inertia in Sees. IVIII Bl and IVIII Dl 
respectively. Sec. VIII C contains a linking discussion. 

are (we believe) novel. In con- 



Sec. 

The results in Sec 



VIII D 



trast, our aim in Sees. VIII A and |VIIIB| is simply to 



reproduce behaviour seen earlier in the literature by lat- 
tice Boltzmann techniques IHl 112] • The purpose of 
this comparative part of our study is threefold: First, and 
foremost, to develop confidence in our own algorithm and 
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the code via which it is implemented; second, to demon- 
strate the method used here, which is potentially simpler 
and faster than lattice Boltzmann, to be equally capa- 
ble of capturing the physics of demixing, in 2D at least; 
and third, to provide an independent check of some re- 
cent results concerning nonequilibrium steady states un- 
der shear |38j. 

In each run, we take as an initial condition at each grid 
point a value of cj) chosen at random from a flat probabil- 
ity distribution between —0.01 and -1-0.01. We checked 
for independence with respect to this initial condition 
(not shown) the statistical properties of the domain mor- 
phologies presented below. This independence holds as 
long as many domains are present, making the system 
self-averaging. Accordingly, for each set of parameters 
we give results below for a single simulation run only. 



A. Zero shear 

Following a temperature quench into the two phase 
regime, domains of each phase form. These progressively 
coarsen in size through the action of surface tension, pass- 
ing through growth regimes that are successively domi- 
nated by diffusive, viscous and inertial dynamics. Within 
these regimes, dimensional analysis predicts a functional 
dependence of the characteristic domain size upon time 
of Ly) cx t^/^, Ly (X t and Li cx t^/^ respectively, as dis- 
cussed in Sec.|VT] Our aim in this section is to investigate 
this behaviour numerically, as done previously by lattice 
Boltzmann methods 32] • 

As noted in Sec. |VII| only a small segment of the full 
coarsening curve L{t) can be explored in any simulation 
run. Accordingly, in different runs we vary the viscosity 
rj and density p. In this way we correspondingly vary 
the crossover length scales Ldv and Lyi. We then com- 
bine the data sets from the different runs into composite 
scaling curves of L/Ldv f^nd L/Lyi, each spanning sev- 
eral decades in scaled length and time. We choose as our 
measure of L the one defined in Eqn.[36|via the structure 
factor. 

To explore the diffusive and viscous hydrodynamic 
regimes, and the crossover between them, we performed 
a single simulation run for each set of parameter values 
of tabl e |(b)[ The composite curve of L/Ldv is shown 
in Fig. l2k). For each segment of the curve, three ma- 
nipulations were performed, (i) The first t = — 40 
time units were discarded to allow for transient dynam- 
ics during the initial process of domain formation, (ii) 
A time-offset tint was subtracted from the time t to al- 
low for these non-universal dynamics during the initial 
transient. In each case this subtraction was performed 
by eye. The appropriate <int was chosen as that which 
gives the most convincing straight line on a log- log plot, 
without any known power present to bias the eye during 
this process, (iii) The data set was cutoff at long times 
once the typical domain size approaches the system size. 
For this measure of the domain size, we used L < 0.02 





FIG. 2: a) Scaled domain size versus scaled time for coars- 
ening in zero shear in the diffusive and viscous hydrodynamic 
regimes. Segments from left to right correspond to runs DVlu 
... DV6u. The vertical arrows show the values of scaled re- 
ciprocal shear rate used in each corresponding run DVls ... 
DV6s. Circles show the times of the top two snapshots of 
Fig. [3] below. 

b) Scaled domain size versus scaled time for coarsening in zero 
shear in the viscous and inertial hydrodynamic regimes. Seg- 
ments from left to right correspond to runs R028u, R022u, 
R029u, R020u, R030u, R019u, R032u. The vertical arrows 
show the values of scaled reciprocal shear rate used in each 
corresponding run R028s ... R032s. Circles show the times of 
the bottom two snapshots of Fig. |3] below. 



as a very conservative criterion ensuring this. The same 
manipulations were performed for runs in table (a) to be 
described below. 

For small values of L/LdVj in the diffusive regime, we 
recover the power Ld oc t^^^ predicted by dimensional 
analysis. This was confirmed previously by lattice Boltz- 
mann simulations J4T]. For large values of L/Ldv, in the 
viscous hydrodynamic regime, dimensional analysis pre- 
dicts the growth law Ly cx t. We instead find L oc t^^^ for 
each run. This anomalous scaling been observed previ- 
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FIG. 3: Snapshots of domain morphology during coarsening 
in zero shear at the times shown by circles in Fig. [2] Top left: 
diffusive regime run DVlu &t t — 1052. Top right: viscous 
hydrodynamic regime run DV5u at t = 131. Bottom left: 
viscous hydrodynamic regime run R022 a,t t — 78.9. Bottom 
right: inertial hydrodynamic regime run R019 at t = 78.9. 



ously, by lattice Boltzmann simulations. It is believed to 
stem from a breakdown of scale invariance in this regime 
in 2D, due to the formation of disconnected droplets [4T] . 
As a result of this breakdown, different measures of L de- 
pend differently on time (not shown here). In 3D (not 
studied here) scale invariance is recovered, along with the 
2/3 growth exponent |42| . 

To study the inertial hydrodynamic regime and (again) 
the viscous hydrodynamic regime, and the crossover be- 
tween them, we use the parameter values of table |(a)[ 
The composite curve of L/Lyi is shown in Fig. [2j3). For 
large values of L/Lyii in the inertial regime, we clearly 
recover the predicted power Li cx t^^^. This has been 
seen before, in three dimensions, by lattice Boltzmann 
simulations [H]. For small values of L/Lyi, in the vis- 
cous hydrodynamic regime, we again find the anomalous 
power L cx t^^^ for each run. This is consistent with 
the breakdown of scale invariance in 2D, and its atten- 
dant departure from the predictions of dimensional anal- 
ysis j41j. Because of this breakdown, segments R028, 
R022, R029 and R029 are not colinear: we do not have 
a single composite scaling curve in this regime. The fact 
that each segment follows a t^^^ power, apparently the 
same as in the inertial regime, is a coincidence resulting 
from having chosen the particular measure of domain size 
given by Eqn. [36] As already noted, the lack of scaling 
in this regime results in different power laws for different 
measures of the domain size 1411. 



sustained over six decades. The same scal- 
7"^/"^ and Lj_ ~ 7^3/4 -^gj-g suggested for 



B. Under shear, with inertia 

We now turn to phase separation in the presence of an 
applied shear flow. The main question of interest here is 
whether shear arrests coarsening and restores a nonequi- 
librium steady state with a finite domain size set by the 
reciprocal shear rate, independent of the system size; or 
whether coarsening persists indefinitely, up to the system 
size, even under shear. 

Despite previous experimental [H [3 [HI |71 HI 13 EHl 
ini Ea Ha [H], numerical [15l [TBI [17l UHl usi EQI EH 
HHIiaillllilllSlllZlllHllMllSnjin] and theoreti- 
cal [T^linjIinilllJIMllMllMJIMllMllMllST] work, this 
question remained open until the recent simulation study 
of Stansell et al. in Ref. [38]. Using lattice Boltzmann 
techniques, they found nonequilibrium steady states of 
the type reproduced by our own simulations in Fig. [6] 
(discussed below). As expected under shear, the domain 
morphology is anisotropic. Accordingly, two length scales 
are needed to characterise it. The domain lengths 
and Ly for the steady states were extracted in Ref. [35] 
via the curvature tensor of Eqn. 37 and scaling plots 
of yj/ Lyi versus l/jTyi were constructed. These 
suggested apparent scaling exponents L^: ~ ^"^/s ^^^^ 
Ly ~ 7" 
ings L|| 

the major and minor principal lengths of the domains. 
Slightly different scalings, discussed below, were obtained 
by the same group in a more recent 3D study |39j . 

In this section we aim to show that our simulations, 
which are in 2D throughout, reproduce these results to 
good approximation. We will thereby gain confidence in 
our technique, before proceeding to the novel contribu- 
tion of this work in Sec. |VIIID| below. Accordingly, we 
perform a single simulation run for each set of parameter 
values used in Ref. [38 , converted into our units. See 
R028, R022, R029, R020, R030, R019 and R032 in ta- 
ble (c) The small discrepancy in the imposed values of 
I/7TV1, evident in Fig.Js] stems from a slightly different 
definition adopted for the interfacial width, realised by 
this author only at a late stage of the present study. 

A typical run took one to two weeks of wall-clock time 
on a Linux box, given exclusive use of a single 3.4GHz 
Intel Xeon CPU with 2Mb cache and 400MHz DDR2 
memory. Approximately 70% of the runtime appears to 
be used switching back and forth between real and recip- 
rocal space at each timestep. 

In each run, we monitored as a function of time the 
characteristic domain sizes Ly and Lj_ extracted from 
the tensor of Eqn. [37] For each of R020, R030, R019 
and R032 we found Ly and Lj_ to saturate at long times, 
showing temporal fluctuations about constant mean val- 
ues (Fig. [4]). For the remainder of this section we focus 
only on these statistically steady states, neglecting the 
first jt — 100 strain units of each run. (This cutoff was 
chosen by a simple visual inspection of Fig. ]4]) In each 
case, finite size effects appear under control, as seen in 
the order parameter snapshots of Fig. ]6] The mean val- 
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0.06 



L 0.04 



0.02 




FIG. 4: Top: larger domain length vs. strain yt. Data sets for 
decreasing values of long term temporal average correspond to 
R020s, R030s, R028b, R029bs, R022bs, R019s, R032s. Bot- 
tom: smaller domain length vs. strain. Data sets for decreas- 
ing values of long term temporal average correspond to RQ20s, 
R030s, R019s, R029bs, RQ32s, R022bs, R028b. Colour online. 



lies of the time series are shown on the scahng plot of 
L/Lyi versus I/7TV1 in Fig. [5] with error bars showing 
the standard deviation. A snapshot of the order param- 
eter for run R032 is shown in Fig. [6] (bottom). 

For runs R028, R022 and R029, in contrast, we found 
the domains to wrap completely round the system in the 
flow direction, eventually comprising trivial horizontal 
stripes connected at the edges of the cell by the periodic 
boundary conditions. We believe this to be due to the 
horizontal system size A^^ being dangerously close to the 
expected values of L|| for these runs, leading eventually 
to nucleation of these stripes. 

To eliminate this effect, we performed new runs R028b, 
R022b and R029b at the same prescribed values of 
I/7TV1 as for R028, R022 and R029, but now for larger 
values of the scaled system size Ax/L^i. This was 
achieved by adjusting Lyi at fixed Aa;,7Tvi. Effectively, 
the left three sets of crosses in Fig. |5]have been slightly 
shifted upward with respect to those in Ref. [38J These 




FIG. 5: Dimensionless scaling plot of lengths vs. shear rate. 
Solid symbols: average of time series of Fig. |4] for L\\,L± for 
strains 'yt > 100. Symbols from left to right correspond to 
R028bs, R022bs, R029bs, R020s, R030s, R019s, R032s. Er- 
rors bars show the standard deviation. Solid lines: power 
law fits to the solid symbols, suggesting Ly ~ ^-O-^is ^nd 
L± ~ ^-0.693 Upper set of crosses shows the scaled sys- 
tem size Ay/Lyi; lower set shows the scaled interfacial width 
l/Lyi. Open symbols: data from Ref. 138! for L\\,L±, repro- 
duced here for comparison by kind permission of the authors 
of that study. 



new runs produced nonequilibrium steady states, as seen 
in the time series of Fig. |4] The long term time averages 
of these series are shown in the scaling plot of Fig. |5] A 
snapshot of the order parameter for run R022b is shown 
in Fig. [6] (top). 

Of course we cannot rule out eventual nucleation of 
system-wrapped stripes at times exceeding those ac- 
cessed numerically. The same comment applies to runs 
R028, R029 and R022 of Ref. [38,. Conversely, the 
system-wrapped stripes seen in the corresponding runs 
in the recent 3D lattice Boltzmann study [H] might well 
be eliminated by adjusting the scaled system size as done 
here. 

Power law fits to our data in the scaling plot 
of Fig. suggest exponents L|| ^ ^-o.6i9±o.oi g^^^^ 
Lj_ r^-o.ra3±o.oo7^ rpj^g quoted uncertainties are the stan- 
dard deviations given by the automated regression pack- 
age, and do not include systematic errors in the simula- 
tions. In contrast Ref. |38j found exponents for Ly , L± 
of -0.678 ± 0.039 and -0.756 ± 0.03, by fitting to the 
open symbols that are reproduced by kind permission in 
Fig. [5] for comparison with our data. It also quoted ex- 
ponents for L^, Ly of -0.678 ±0.042 and -0.759 ±0.029. 
The more recent 3D study |39j by the same group found 
exponents for Ly, L_l of -0.64 ± 0.06 and -0.67 ± 0.03, 
with —0.53 ± 0.04 in the third dimension. Within their 
0(10%) margins of error, these agree with the exponents 
for Lu and Lj_ found in the 2D study of the present work. 
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FIG. 6: Snapshots of the steady state order parameter for 
R022b at strain 'jt = 112 (top) and R032 at strain 'jt = 161 
(bottom). 



C. Discussion 

In this section, wc will discuss in more detail the shape 
of the scaling plot of L/Lyi versus If-yTyi in Fig. [s] In 
doing so, we will motivate the further numerical study of 
Sec. IVniDJ below. 

Translated into equivalent scaled times t/Tyi, the 
range of scaled shear rates If-^Tyi explored in Fig. [s] 
would span the viscous and inertial hydrodynamic 
regimes in the coarsening plot of Fig.l2b). One possibility 
for the shape of the scaling curve L/Lyi versus 1/jTyi 
under shear is that it should be the same as that of the 
corresponding coarsening plot of L/Lyi versus i/Tvi in 
zero shear, given the simple substitution t ~ 7^^. In view 
of this, it is instructive to compare Fig. [sjwith Fig. ^p) 
in some detail. 

In the coarsening case (Fig. [2]3), two distinct regimes 
are apparent. The inertial hydrodynamic regime, on the 
right side, shows the t'^^^ exponent predicted by dimen- 
sional analysis. Each segment neatly lines up with the 
next, consistent with dynamical scaling. In contrast, in 
the viscous hydrodynamic regime on the left hand side 
the segments do not align. Furthermore, each departs 
from the predicted power. As noted above, these 
anomalous features stem from non-scaling effects that 
arise in 2D. In 3D these are suppressed and the pre- 
dicted power law is recovered. The counterpart of 
Fig. [2|d) for 3D systems thus comprises a clean in the 
viscous regime on the left hand side, and t^^^ in the in- 
ertial regime on the right hand side, with a rather slow 
crossover in between. See Fig. 9 of Ref. |12]. An ap- 



plied shear flow is also anticipated to suppress these non- 
scaling effects. Accordingly, the 2D results of Fig. |5]are 
expected capture the basic physics and, indeed, compare 
favourably with the recent 3D sheared study of Ref. [12] . 

With these remarks in mind, it is now clear that we 
should in fact compare Fig. [5] for sheared systems with 
Fig. 9 of Ref. [42 for coarsening systems. As noted 
above, one might naively expect the two plots to have 
the same shape, up to the substitution t = 7^^. 

Instead, two differences are clearly apparent. First, 
Fig. |5]has two characteristic lengths, in contrast to the 
single length that characterises unsheared systems in 
Fig. ^p). With hindsight this is obvious: the domain 
morphology is anisotropic in sheared systems, and so 
characterised by two lengths. 

The second, and more subtle, difference is that sheared 
systems apparently lack a distinction between viscous 
and inertial regimes, with a single power law spanning 
all six decades in Fig. |5] A possible explanation for this 
is that the range of scaled shear rates in Fig. [5] actually 
occupies a window of extremely slow crossover between 
a true viscous regime at smaller I/7TV1 and a true in- 
ertial regime at larger 1/jTyi, reminiscent of the slow 
crossover seen in the coarsening plot in Ref. [42 . An- 
other is that the value of l/7rvi at which the crossover 
occurs in the sheared case is much smaller than the corre- 
sponding crossover value of t/Tyi in the unsheared case, 
which is formally 0(1) but in practice lies around lO''. 
The results in Fig. [s] would then all lie in the inertial 
regime, to the right of this crossover. 

To test these ideas, in the next section we perform 
simulations at strictly zero Reynolds number, p = 0. 
We thereby take the limit 7TV1 — > 00 at the outset and 
perform simulations for several finite values of 1/7Tdv- 
When converted into equivalent scaled times t/Toy these 
span the diffusive and viscous hydrodynamic regimes, as 
shown by the vertical arrows in Fig.|2^). By analogy with 
the fact that the 3D counterpart of the viscous regime of 
Fig. |2^) would be expected in the limit i/Tov ^ 00 to 
match into that of the viscous regime in Fig. 9 of Ref. [15] 
in the limit i/TVi — + 0, we might then expect to access, 
for large I/7TDV ~* 00, the scalings that Ly and L± 
would attain in the limit 1 /"fTyi — )■ at the extreme left 
of an enlarged Fig. |5] 

Instead, however, we will find no evidence for non 
equilibrium steady states in inertialess systems. In con- 
trast, coarsening appears to persist indefinitely, up to 
the system size, despite the applied shear flow. Given 
the existence of non equilibrium steady states for values 
of I/7TV1 at the left hand edge of Fig. [s] where the ef- 
fects of inertia are non-zero but likely to be small, this 
suggests that inertia plays the role of a singular pertur- 
bation in this problem. Indeed, we will argue in Sec. |IX| 
below that the single scaling seen for each of L|| and L± 
across Fig. [s] results from a mixed visco-inertial regime 
across the entire plot. The crisis of infinite aspect ratio 
L\\/L±_ — > CO in the hmit l/7/Tvi 0, suggested by 
the data in Fig. [5] is consistent with our suggestion that 
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FIG. 7: Top: larger domain length L\\ vs. strain -yt. Bot- 
tom: smaller domain length L± vs. strain jt. In each case, 
decreasing values of L at fixed 7t = 15 correspond to runs 
DV5s, DV6s, DV4s, DV3s, DV2s, DVls. 

inertia plays a singular role. 

D. Under shear, without inertia 

Motivated by the discussion of the previous section we 
now consider phase separation under shear in inertialess 
systems, setting p = 0. As just discussed, by studying 
the limit t/Tiyy — > cx) we might then, a priori, have ex- 
pected to gain some insight into the far left hand side of 
an enlarged Fig. [s] l/7rvi ^ 0. Accordingly, we per- 
form a single simulation run for each of the parameter 
sets DVls to DV6s in table |(d)[ The values of the scaled 
reciprocal shear rates 1 /'jTdv are marked as correspond- 
ing scaled times </Idv on the coarsening plot of Fig. [2^), 
where they are seen to span both the diffusive and viscous 
hydrodynamic regimes. 

In each run we monitored the characteristic domain 
sizes L|| and L± as a function of time. As can be seen in 
Fig. [TI in each case the larger length Ly grows without 
bouiia until it attains the system size 0(1). Correspond- 
ingly, a snapshot of the order parameter after many strain 
units reveals system- wrapped domains, with pronounced 
finite-size effects. See Fig. [8] The snapshots for DVls 



to DV3s strongly resemble those reported in Ref. (3T] for 
model B, which lacks hydrodynamics. This is consistent 
with the fact that these runs occupy the diffusive regime 
when marked as equivalent scaled times in Fig.|2|D). 

Beyond the parameter sets of table |(d)[ we fur- 
thermore performed runs for shear rates 7 = 
0.1,0.03,0.01,0.003 at each value of the viscosity rj — 
1000.0,100.0,10.0,1.0,0.3,0.1, thereby covering a com- 
plete rectangle in (ry, 7)-space, in contrast to the single 
slice taken by DVls to DV6s. (For historical reasons 
these runs had a slightly larger value I = 0.0025 for the 
interfacial width, but we do not expect this to make a 
qualitative difference.) We found no evidence in any run 
for a nonequilibrium steady state, unlimited by finite-size 
effects. 

These results clearly suggest that the limit I/7TDV 
00, which is approximated by runs DV5s and DV6s, does 
not match the limit I/7TV1 — > 0, which is approximated 
by runs R028b, R022b. In the next section, we propose 
that this is because the limit If-^TDy 00 corresponds 
to a pure viscous regime in which no steady state exists; 
while the limit 1 / jTyi — > corresponds to a mixed visco- 
inertial regime in which a steady state does exist. 



IX. ABSENCE OF NONEQUILIBRIUM 
STEADY STATES IN INERTIALESS SYSTEMS 

We now present an analytical argument in support of 
the above numerical observation: that nonequilibrium 
steady states are not attained in inertialess systems. As 
we shall find, a closely related argument supports the 
existence of a single power law for each of Ly , L± across 
the whole of Fig. [s] as seen first in Ref. [SS]. Recall that, 
a priori, we might have expected this plot to comprise 
separate viscous and inertial regimes. 

As a preliminary step, we recall the case of domain 
coarsening in an unsheared system. As discussed above, 
as time proceeds the system passes through regimes that 
are successively dominated by diffusive, viscous and in- 
ertial dynamics. In each regime, the system is aware of 
the surface tension cr; the time t; and either the mobility 
M (diffusive regime), the viscosity 77 (viscous regime) or 
the density p (inertial regime). In each case, there exists 
only one possible combination of these relevant parame- 
ters that has the dimensions of a length. Thus we have 
the following predicted growth laws for the characteristic 
domain size: 

t / t"^ \ ^^"^ 

Lu = {Matf^, Lv = — , and Li = ( — I . (39) 
V \ P J 

Consider next a nonequilibrium steady state under 
shear. Neglecting fiuctuations, the time t is now irrel- 
evant. In its place, however, we gain the reciprocal shear 
rate as a relevant parameter. If "pure" diffusive, viscous 
and inertial regimes were still to exist, then we could 
again only construct a single length scale to characterise 
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(d)DV4s, 7i = 63.2 (e)DV5s, -yt = 55.3 (f)DV6s, -yt = 47.4 



FIG. 8: Snapshot domain morphologies after many strain units 3 = ^1 for reciprocal shear rates that, as equivalent times 7 
would be shown by the arrows in Fig. [2k). These runs have zero inertia, p = 0. 



each: 



(77 



and Lts = 



(77 



1/3 



(40) 

by direct comparison with (39). 

In contrast, however, our numerical results clearly 
show the domain morphology to be anisotropic, as ex- 
pected in a sheared system. It is therefore characterised 
by two lengths, which (according to Fig. [s]) scale differ- 
ently with the applied shear rate. A steady state under 
shear therefore cannot exist in a "pure" diffusive, viscous 
or inertial regime, because in each there are insufficient 
parameters out of which to construct the two lengthscales 
needed to characterise it. 

How can we retain enough parameters, out of the can- 
didates cr, 7, M, p, T] and t, to construct the two length- 
scales needed to characterise an anisotropic domain mor- 
phology under shear? Assuming that cr and 7 are always 
relevant, two options are as follows. 



1. If the system is to exist in a pure diffusive, vis- 
cous or inertial regime, with knowledge of just one 
of M, 7y or p (as well as a and 7), it must retain 
dependence on time t. Therefore, it must fail to 
attain a steady state. 

Assuming power laws, a purely diffusive regime 
would then have 



Lr 



{Maty/%itr, 



(41) 



in which a assumes two different values, which 
would be prescribed by more detailed physics than 
is contained in the present argument. Likewise, a 
purely viscous regime would have 

Lvs = -m\ (42) 



again with two different values for b. 

2. If the system is to attain a steady state, thereby 
losing knowledge of the time t, it must retain de- 
pendence upon at least two of M, rj and p. For 
example, a shear-induced steady state that is free 
of diffusion on the lengthscale of domains (no M) 
must exist in a mixed viscous-inertial regime with 
knowledge of both viscosity 77 and density p: 



L 



vis 



(43) 



with two different values for c. 



We propose that both of these cases are seen in our 
numerical simulations. When inertia is present, for finite 
values of p (however small) , the system exists in a mixed 
visco-inertial steady state, as discussed in option 2. This 
is consistent with the observation of a single power law 
scaling for each of iy and L± across all six decades in the 
plot of L/Lyi versus I/7TV1 in Fig. [s] and in Ref. 38J. In 
contrast, when inertia is strictly absent (p = 0, infinite 
7rvi)i the system exists in a pure diffusive or viscous 
regime and never attains a steady state, as in option 1. 

To summarise: we suggest that the limit I/t'Tdv ^ 
00 corresponds to a pure viscous regime with no steady 
state, while the limit I/7TV1 — *■ corresponds to a mixed 
visco-inertial steady state. Thus we propose finally that 
inertia provides the role of a singular perturbation in this 
problem. 



X. SUMMARY AND OUTLOOK 

We have studied numerically phase separation in bi- 
nary fluids with full hydrodynamics in two dimensions. 
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considering both (1) unsheared and (2) sheared systems, 
both (a) with and (b) without inertia. Of these, cases 
(la), (lb) and (2a) have been studied by previous authors 
using Lattice Boltzmann (LB) methods pgl 155 1 liTl li^ . 
In this paper, we have introduced an alternative simu- 
lation technique that uses finite-differencing and spec- 
tral methods. We have also used a convenient trans- 
formation to render trivial the implementation of Lees- 
Edwards sheared periodic boundary conditions jiO] . 

In unsheared systems, phase separation occurs via a 
process of domain coarsening. Our simulation method 
successfully recovers results obtained previously by LB 
for this process, both with and without inertia (la, lb 
above). In particular, it finds the famihar power law 
Lq ~ t^^^ characterising the growth of the typical do- 
main size in the diffusive regime |41j . It also recovers 
L\ ^ t^l'^ in the inertial hydrodynamic regime [42]. In 
the viscous hydrodynamic regime it finds the anoma- 
lous power L ^ t^/'^, compared to the predicted one of 
iv ^ t^- As noted by previous authors, this is due to 
subtle non scaling effects that arise in 2D from the for- 
mation of disconnected droplets, also seen in LB stud- 
ies IH]. Such effects are eliminated in 3D [15] and also 
seem suppressed under shear [251 1^ . 

We have also successfully reproduced the observations 
of existing LB studies for sheared systems that have non- 
zero inertia (case 2a above) |38j. Here, an applied shear 
flow arrests domain coarsening and restores a noncqui- 
librium steady state with domains of a finite size set 
by the inverse shear rate. The domain morphology is 
anisotropic, characterised by two lengthscales 
Scaling exponents Ly ~ ^-o.eig g^^^ ^ ^-0.693 fQ^^d 
here agree with those of LB, to within margins of error. 
An outstanding puzzle, however, is why these 2D expo- 
nents agree better with the exponents Ly ~ ^-0-64^ ^ 
7-°-^^ of the 3D LB study [39], than those of the 2D LB 
study [35;, i|| - 7-0-678 j^jj^ _ ^-0.756^ rj,^ investi- 
gate this, it would be interesting to study the role of sys- 
tematic errors in our simulations; to consider the possible 
eventual nucleation of system wrapping stripes in the 2D 
LB study of Ref. [38 for the smaller values of I/7TV1; 
and to extend the present finite-differencing work to 3D. 

Our successful recovery of these important existing 
results in regimes (la, lb, 2a) provides some confi- 
dence in our simulation method. Beyond thereby hav- 
ing demonstrated this method to be capable of captur- 
ing the physics of demixing, the other main contribu- 
tion of this paper has been a novel study of phase sep- 
aration in sheared systems that are strictly inertialess 
(case 2b) . Here we found no evidence of nonequilibrium 
steady states, free of pronounced finite size effects. In- 
stead coarsening appears to persist indefinitely until the 
typical domain size attains the system size, as in zero 
shear. 

To support this observation, we have suggested by 
means of a simple analytical argument that sheared in- 
ertialess systems adopt either a pure diffusive or pure 
viscous regime, in each of which there are insufficient pa- 



rameters out of which to construct the two lengthscales 
needed to characterise an anisotropic domain morphol- 
ogy in steady state. By extending this argument slightly, 
we have also suggested that sheared systems with any 
amount of inertia, however small, exist in a mixed visco- 
inertial steady state. This provides a possible explana- 
tion for the observation of a single scaling with shear rate 
for each of and L± across all six decades in Fig.jSj and 
in the corresponding plot of Ref. [35] 

If this suggestion is correct, it remains unclear why the 
viscous and inertial regimes should mix to yield a steady 
state, while the diffusive and viscous regimes apparently 
remain separate, precluding nonequilibrium steady states 
in truly inertialess systems. A possible explanation lies in 
the more severe nonlinearity (in v) of the inertial terms 
in the equation of motion. 

In future work, we aim to investigate whether the ab- 
sence of nonequilibrium steady states in inertialess sys- 
tems persists in 3D. In extending our method to 3D, 
several challenges are to be faced. Of these, the main 
ones appear to be contained in the basic Navier-Stokes 
equations of incompressible fluid flow, and are not com- 
plicated significantly by any additional order parame- 
ters [48]: (p in this model. However, relatively stan- 
dard methods do exist for finite-differencing the incom- 
pressible Navier-Stokes equations in 3D, as discussed 
in Ref. [iS]. At each timestep these involve updating 
the vorticity via a slightly modified equation of motion 
dtOJ — ■ ■ ■ then updating the velocity field either in the 
velocity- vorticity formulation (as effectively done here in 
2D via the intermediate of the streamfunction) or in the 
vector potential-vorticity formulation. In neither case 
does the pressure field need to be calculated directly. An 
outstanding issue before any such algorithm could be run 
efficiently concerns its scaling with system size, and its 
corresponding ease of parallelisation. The same question 
is also relevant in 2D, and a direct comparison of our 
algorithm with that of Ref. [38j would clearly be inter- 
esting. In the longer term, the outcome of such stud- 
ies might help determine whether finite differencing can 
emerge as a useful tool in such problems, alongside the 
already tested and reliable LB methods. 

Other open questions concern the role of initial condi- 
tions in sheared systems. All the simulations reported 
here consider a temperature quench performed in the 
presence of a shear flow. Future work should consider 
an already demixed system, with either a flat interface 
or a minimal-surface droplet, subsequently subject to a 
sudden shear startup. We also aim to address demixing 
in complex macromolecular fluids, focusing on the role of 
viscoelasticity. 

Acknowledgements 

The author thanks Prof. Alan Bray, Prof. Mike Gates 
and Prof. Peter SoUich for helpful discussions and feed- 
back, and the UK's EPSRC GR/S29560/02 for funding. 



14 



APPENDIX A: TRANSFORMATION TO THE COSHEARED FRAME 



Here we give details of the transformation to the cosheared frame [4D, . As discussed in Sec. |v] above, this is 
performed in order to render trivial the numerical implementation of sheared periodic boundary conditions. As a first 
step, we separate the velocity field into a constant affine contribution '^yx and a fluctuating part v 



v{x, t) — jyx + v{x, t). 



(Al) 



Noting that jyx automatically satisfies the incompressibility condition, we define the fluctuating parts of the stream- 
function and vorticity via 



and 



■u = V A tpz, 



Our equation set then comprises ( A2 ) and ( A3 1 together with 

p{dtOj + ■u.V ij) + pjy dxi^ — vj^I^uj — [V A 0V^] • i, 

from Eqn. |3] 



a* + i3. V + 72/ d^ct) = ?2v2^^ 



from Eqn. |7] and 



1)-PV 



2v72, 



from Eqn. |8] unchanged. 

We then make a transformation to the cosheared frame 

(x, y, t) ^ {x ^x- jty, y = y, t' = t). 

The various partial derivatives then become 

{dx,dy, dt) = {dx',-jtdx' + dy',-jydx' + df). 

Accordingly, we define the cosheared 2D gradient operator 

Vc = xdx' +y{-itdx' +dy'). 

Finally, for any function a we write 

a{x,y,t) = A{x',y',t'). 



(A2) 
(A3) 
(A4) 
(A5) 
(A6) 

(A7) 
(A8) 
(A9) 
(AlO) 



Throughout we continue to work with velocity components , Vy S ind not Vxi,Vyi. 

Inserting Eqn. A2 into the v ■ V terms on the LHS of Eqns. A4 and A5 and performing the transformation (A7l 
to (AlO I on Eqns. A3 to A6 we get the equation set 



together with 



pdt'n + p {dy'^d^.n - d^'^dy'n^ = V^f^ - [dx:<^dy-M - dy.^d^-M] 



dt'<^+[dy:^dx^<^-dx.^dy.<^] ^fvlM, 



(All) 

(A12) 
(A13) 

(A14) 

In these M represents upper case p, not mobility. A priori, the bracketed expressions in Eqns. A12 and A13| contain 
terms in the applied shear rate 7. However in each bracket these are equal and opposite, and so cancel. For clarity 
we finally drop the tildes and dashes, and revert from upper to lower case. The final governing equations are then as 
summarised in Sec. |V] above. 



( 



and 



M = $($^ - 1) - rV^$ 
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APPENDIX B: NUMERICAL METHOD 



Here we discuss the numerical algorithm used to study the dynamical evolution of <j){x,t),Lu{x,t) and ip{x,t), as 
specified by Eqns. [21] to [24] and [26] Our basic strategy is to step along a grid of time values = nAt for n = 1, 2, 3 • • • . 
Discretization with respect to time of any quantity / is denoted /(t") = or sometimes by /|". At each timestep, we 

in three separate stages. First, we update the compositional order parameter 



update 



n+l „l,n+l 



n+1 



according to Eqns. 23 



and 



24 



with fixed, old values of the stream- function ip^. We then update 
using Eqn. 22 at fixed (j)"^^, ■(/'"• Finally we update the streamfunction -0" 7/)"+^ using Eqn. 21 at fixed lu^' 



,n+l 



1. The update 1/)"+^ using Eqns. 23 and 24 is performed in two successive partial updates. In the first we 

implement the advective term in Eqn. 23 to give 0" 0"+^/^. In the second we implement the diffusive term to 
j^n+i _ fpj-^g advective term is handled using an explicit Euler algorithm [5(7. Temporarily setting 



give 



in+1/2 



aside the issue of spatial discretization, this can be written 



(Bl) 



This is then spatially discretized on a rectangular grid of (A^N/Ay) x N mesh points in real space (x — y). 



with constant mesh intervals 6x 



Ay/N. Using indices i = I ■ ■ ■ A^N/Ay and j = 1 • • • iV, we denote any 



discretized function f{xi,yj) = fij, or sometimes f\ij. Periodic boundary conditions are imposed by setting 



— fi(N-l), fiO — fiN, fi{N+l) — fill fi(N+2) 

in Eqn. |Bl|are discretized as follows: 



1 

26x 



fi2, and similarly in the x direction. The derivatives of 4' 

(B2) 



with 



and 



25„ { ( 



1 



1 



^^0-1)) + ^ (V'(»-2)0+i) - V'(i+2)0--i))} if s > 0, 



(B3) 



(B4) 



(B5) 



The derivatives of in Eqn. Bl are discretized in the same way. 
It then remains to implement the diffusive part of Eqn. |23| 

with f — (f)^ . After calculating / on our rectangular grid in real space, this equation is handled in reciprocal space 
by taking fast Fourier transforms x ^ and y Qy using a standard NAG routine |51j . The transformation 
in each dimension generates a real and an imaginary part, so for each mode q = (qx, qy) we need to consider a 
vector of the (transposed) form cj)^ — {(f)^^, ^ri, (jiu) where subscript "r" denotes real part, and "i" imaginary. 
The respective transforms D2 and i?4 of the operators P and l'^ can easily be found analytically: 



D2 = V 



( B -A \ 

£» A 

A £> 

V -A D ) 



and Da = I 






V -27:>A 





D^ + A^ 
2DA 




-2L>A \ 
2DA 
D^ + A^ 

D^ + A^/ 



in which 



D^-{aql + ql), A ^ bq^qy with a = 1 + [s(^)]^ b = 2s{t). 
For each q-mode, we then have 

dt4)^-D2-c()-Di-4>+D2- f. 



(B6) 



(B7) 



(B8) 



To evolve this in time, we use an explicit Euler algorithm for the first and third terms, and a semi-implicit 
Crank-Nicolson algorithm for the second term. Thus we have 



in which Dr 



^n+l _ ^n+1/2 ^ . ^n+1/2 _ 1^^ . ^ ^n+1/2^ ^ fj^ . j^n ^ 

St Dm for m = 2, 4. Rearranging gives finally 



^"+1 = 0" 



-1/2 



-1/2 



D2 f 



(B9) 



(BIO) 
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2. We now update uj^~^^ using Eqn. 22 at fixed cj)'^'^^ , . The advective term on the LHS is updated in 

the same way as its counterpart in the (j) equation above. To avoid ineflicient multiple switching between real 
and Fourier space, this is in fact done at the same time as the corresponding update of 4> in 1. above. (This 
reordering leaves the algorithm exactly unchanged.) We then update the RHS of Eqn. 22 Divided across by p, 
this reads: 



in which v ~ rj/p and 



9 



[dx(t)dyp. - dy4>dxfi] 



As a first step, we calculate g{x, y) using the newly updated from 1. To do so, we first calculate 

p = - 1) - f{l + [s{t)]^)dlcb - l^dl^p + 2fsd,dycl). 
At each grid point we spatially discretize the derivatives in this expression according to 

1 



similarly for d'^<f>, and 



dxdy4>\ij 



^5x5y 



(i-l)0-l)J 



(Bll) 



(B12) 



(B13) 



(B14) 



(B15) 



The first order derivatives of and p with respect to x and y in (B12l are then discretized as in (B2|. 
Having calculated g{x,y) in real space, we then Fourier transform Eqn. Bll to get 

dtLj = C ■ uj + g, 



(B16) 



in which we have used the same vector/matrix notation as in 1 above, with C = 1^02/1^. To evolve this in time, 
we use a semi-implicit Crank-Nicolson algorithm for the first term on the RHS to get 



,n+l 



(B17) 



with C — StC and g = Stg. The superscript n + 1 on the last term (B17l serves to remind us that g was 
calculated using the new 0"+^ from part 1. Rearranging, we get finally 



n+l 



w" + (^-iC)-i-(C 



(B18) 



3. We finally update the streamfunction -0" -0"+^ using Eqn. 21 at fixed 0;"+^. For each <7— mode, we have 

(B19) 



n+l 



in which E = 



£)2 



[do a\ 

D -A 

-A £> 

V A D / 



Finally, we transform all functions back to real space and return to step 1 to start the next timestep. 

Algorithm at zero Reynolds number 

The algorithm discussed so far is suited to non-zero values of the fluid density p. At zero Reynolds number, with 
p = 0, Eqn. [21] and |22| of our basic equation set combine to give the simpler equation 



= -?7 VcV' - [dx(j)dyp - dycjidxp] , 



(B20) 



and we need no longer consider the vorticity field uj. Equations |23| and |24| remain unchanged. Correspondingly, step 
1 of our algorithm is also unchanged. Steps 2 and 3 now combine to give 



n+l 



E.g 



n+l 



(B21) 
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As above, we calculate g in real space then take a Fourier transform. In Fourier space 

/ L>2 + A2 2D A \ 

D^ + A^ -~2DA 

-2DA D'^ + A^ 

\ 2D A D^ + A^ J 



1 



+ A2)2 - {2DAy 



(B22) 



After calculating ip, we revert to real space to start the next timestep. 
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